Assessing the impact of free-roaming dog population management through systems modelling

Free-roaming dogs can present significant challenges to public health, wildlife conservation, and livestock production. Free-roaming dogs may also experience poor health and welfare. Dog population management is widely conducted to mitigate these issues. To ensure efficient use of resources, it is critical that effective, cost-efficient, and high-welfare strategies are identified. The dog population comprises distinct subpopulations characterised by their restriction status and level of ownership, but the assessment of dog population management often fails to consider the impact of the interaction between subpopulations on management success. We present a system dynamics model that incorporates an interactive and dynamic system of dog subpopulations. Methods incorporating both fertility control and responsible ownership interventions (leading to a reduction in abandonment and roaming of owned dogs, and an increase in shelter adoptions) have the greatest potential to reduce free-roaming dog population sizes over longer periods of time, whilst being cost-effective and improving overall welfare. We suggest that future management should be applied at high levels of coverage and should target all sources of population increase, such as abandonment, births, and owners of free-roaming dogs, to ensure effective and cost-efficient reduction in free-roaming dog numbers.

There are between 700 million and one billion dogs globally, comprising owned and unowned dogs that live in shelters, homes, or freely in urban and rural environments 1,2 . Where dogs freely roam and exist in high densities, they can present significant challenges to public health 3,4 , wildlife conservation 5 , and livestock production 6,7 . Free-roaming dogs (i.e. unrestricted dogs found on streets) can also experience conditions leading to poor health and welfare states 8,9 . Population management is conducted to control free-roaming dog population size and, depending on the approach taken, to improve dog health and welfare and mitigate public health and conservation problems 10,11 .
Dog populations can be divided into subpopulations depending on their relationship with humans (owned or unowned) and their restriction status (restricted or free-roaming). The free-roaming dog population (including both owned and unowned free-roaming dogs) is of most interest for dog population management 12 . Dog subpopulations are interactive and dynamic-individuals can move between different subpopulations through time (Fig. 1). Throughout its life, a dog could spend time in multiple states-for example, living as an unowned free-roaming dog, being caught and moved into shelter accommodation, before being adopted into a home.
Dog population management can involve reproductive control (e.g. catch-neuter-release; CNR), sheltering, culling, and responsible ownership campaigns 10,12 . The free-roaming dog population size can be reduced by decreasing recruitment into (e.g. births and immigration) and increasing removal from (e.g. mortality and emigration) the population. As the dog population consists of interacting subpopulations, migration not only occurs between free-roaming dog populations in different geographic locations, but also between the subpopulations within a location, providing additional sources of recruitment. In terms of their effect on recruitment and removal, both culling and sheltering aim to increase rates of removal from the free-roaming dog population. Both CNR and responsible ownership campaigns aim to decrease rates of recruitment to the free-roaming dog population: CNR through decreased birth rates; and responsible ownership campaigns, depending on the approach, through the reduced abandonment and roaming of owned dogs.
Modelling methods can aid in planning, implementing, and evaluating dog population management strategies. System dynamics modelling allows complex and interactive systems to be explained and the impact of interventions on model behaviour to be evaluated 13,14 . Dog population dynamics have been investigated with both system dynamics approaches [15][16][17][18][19][20][21][22][23] and agent-based models 24,25 , assessing the effects of management strategies on population size 15,21,22,25,26 , euthanasia rates 23 , and disease dynamics 19,20 . Most previous studies have modelled dynamics within a single subset of the population (e.g. owned or both owned and unowned free-roaming). Few have modelled the interactions between these subpopulations. Modelling these interactions could provide greater insight into the effectiveness of dog population management interventions in reducing dog population size. All population management methods require the investment of resources (e.g. staff, facilities, and equipment), but few previous models have considered cost-effectiveness 16,22,23 .
The aim of this study was to compare the impact of different dog population management methods on (i) the sizes of the different dog subpopulations, (ii) the staff-resources required, and (iii) the welfare costs of each intervention using a system dynamics modelling approach. The management methods under investigation included sheltering, culling, CNR, and responsible ownership. To reflect that management practices often involve combinations of management methods 27 , we also investigated combined CNR and sheltering, and combined CNR and responsible ownership.

Methods
Model description. The system dynamics model divided an urban dog population into the following subpopulations: (i) free-roaming dogs (both owned and unowned free-roaming, i.e. unrestricted dogs found on streets), (ii) shelter dogs (unowned restricted dogs living in shelters), and (iii) owned dogs (owned home-dwelling restricted dogs) (Fig. 1). The subpopulations change in size by individuals flowing between the different subpopulations or from flows extrinsically modelled (i.e. flows from subpopulations not included in the systems model; the acquisition of dogs from breeders and friends to the owned dog population, and the immigration/ emigration of dogs from other neighbourhoods).
Ordinary differential equations were used to describe the dog population dynamics. The models were written in R version 3.6.1 28 , and numerically solved using the Runge-Kutta fourth order integration scheme with a 0.01 step sizes using the package "deSolve" 29,30 . For the baseline model, Eqs. (1-3) were used to describe the rates of change of dog subpopulations in the absence of management.
Baseline free-roaming dog population (S): www.nature.com/scientificreports/ In the baseline model, the free-roaming dog population (Eq. 1) increases through the free-roaming dog intrinsic growth rate (r s ), and the abandonment and roaming of dogs from the owned dog population (α) and decreases through adoption to the owned dog population (δ). The intrinsic growth rate is the sum of the effects of births, deaths, immigration, and emigration, which are not modelled separately. In this model, the growth rate of the free-roaming dog population is reduced depending on the population size in relation to the carrying capacity, through the logistic equation (r real = r max (1 − S/K s )) 31 . In the baseline simulation, the free-roaming dog population rises over time, until it stabilises at an equilibrium size.
Baseline shelter dog population (H): The shelter dog population (Eq. 2) increases through relinquishment of owned dogs (γ) and decreases through the adoption of shelter dogs to the owned dog population (β), and through the shelter dog death rate (µ h ). There is no carrying capacity for the shelter dog population as we assumed that more housing would be created as the population increases. This allowed calculation of the resources required to house shelter dogs.
Baseline owned dog population (O), The owned dog population (Eq. 3) increases through the owned dog growth rate (r o ), adoption of shelter dogs (β), and adoption of free-roaming dogs (δ); and decreases through abandonment/roaming (α) and relinquishment (γ) of owned dogs to the shelter dog population. The growth rate of the owned dog population (r o ) combines the birth, death, and acquisition rates from sources other than the street or shelters (e.g. breeders, friends) and was modelled as density dependent by the limit to growth logistic formula (1 − O/K o ).
Parameter estimates. Detailed descriptions of parameter estimates are provided in the supplementary information. The simulated environment was based on the city of Lviv, Ukraine. This city has an area of 182 km 2 and a human population size of 717,803. Parameters were estimated from literature, where possible, and converted to monthly rates (Table 1). Initial sizes of the dog populations were estimated for the baseline simulation, based on our previous research in Lviv 32 . The carrying capacity depends on the availability of resources (i.e. food, shelter, water, and human attitudes and behaviour 33 ) and is challenging to estimate. We assumed the initial freeroaming and owned dog populations were at carrying capacity. Initial population sizes for simulations including interventions were determined by the equilibrium population sizes from the baseline simulation (i.e. the stable population size, the points at which the populations were no longer increasing/decreasing).
Estimating the rate at which owned dogs are abandoned is difficult, as abandonment rates are often reported per dog-owning lifetime 32,34 and owners are likely to under-report abandonment of dogs. Similarly, it is www.nature.com/scientificreports/ challenging to estimate the rate that owned dogs move from restricted to unrestricted (i.e. free-roaming). For simplicity, we modelled a combined abandonment/roaming rate (α) of 0.003 per month, estimated based on our previous research in Lviv and from literature [34][35][36] . We derive the owned dog relinquishment rate (γ) from New et al. 37 . We estimated shelter (β) and free-roaming adoption rates (δ) from shelter data in Lviv. We set the maximum intrinsic growth rate for the free-roaming dogs (r s ) at 0.03 per month, similar to that reported in literature 17,19,38 . We assumed that demand for dogs was met quickly through a supply of dogs from births, breeders and friends and set a higher growth rate for the owned dog population (r o ) at 0.07 per month. We assumed shelters operated with a "no-kill" policy (i.e. dogs were not killed in shelters as part of population management) and included a shelter dog death rate (µ h ) of 0.008 per month to incorporate deaths due to euthanasia for behavioural problems or health problems, or natural mortality. We modelled neutered free-roaming dog death rate (µ n ) explicitly for the CNR intervention at a minimum death rate of 0.02 per month 38-41 . Interventions. Six intervention scenarios were modelled ( Table 2): sheltering; culling; CNR; responsible ownership; combined CNR and responsible ownership; and combined CNR and sheltering, representing interventions feasibly applied and often reported 27 . Table 2 outlines the equations describing each intervention. To simulate a sheltering intervention, a proportion of the free-roaming dog population was removed and added to the shelter dog population at sheltering rate (σ). In culling interventions, a proportion of the free-roaming dog population was removed through culling (χ).
Free-roaming dog population with sheltering intervention: Shelter dog population with sheltering intervention: Free-roaming dog population with a culling intervention: To simulate a CNR intervention, an additional subpopulation was added to the system (Eq. 7): (iv) the neutered free-roaming dog population (N; neutered, free-roaming). In this simulation, a proportion of the intact (I) free-roaming dog population was removed and added to the neutered free-roaming dog population. A neutering rate (φ) was added to the differential equations describing the intact free-roaming and the neutered free-roaming dog populations. Neutering was assumed to be lifelong (e.g. gonadectomy); a neutered free-roaming dog could not re-enter the intact free-roaming dog subpopulation. Neutered free-roaming dogs were removed from the population through the density dependent neutered dog death rate (µ n ); death rate increased when the population was closer to the carrying capacity. The death rate was a non-linear function of population size www.nature.com/scientificreports/ and carrying capacity modelled using a table lookup function (Fig. S1). Neutered free-roaming dogs were also removed through adoption to the owned dog population, and we assumed that adoption rates did not vary between neutered and intact free-roaming dogs. Neutered free-roaming dog population: Intact free-roaming dog population with neutering intervention.
To simulate a responsible ownership intervention, the baseline model was applied with decreased rate of abandonment/roaming (α) and increased rate of shelter adoption (β). To simulate combined CNR and responsible ownership, a proportion of the intact free-roaming dog population was removed through the neutering rate (φ), abandonments/roaming decreased (α) and shelter adoptions increased (β). In combined CNR and sheltering interventions, a proportion of the intact free-roaming dog population (I) was removed through neutering (φ) and added to the neutered free-roaming dog population (N), and a proportion was removed through sheltering (σ) and added to the shelter dog population (H).
Intact free-roaming dog population with combined CNR and sheltering interventions: Intervention length, periodicity, and coverage. All simulations were run for 70 years to allow populations to reach equilibrium. It is important to note that this is a theoretical model; running the simulations for 70 years allows us to compare the interventions, but does not accurately predict the size of the dog subpopulations over this long time period. Interventions were applied for two lengths of time: (i) the full 70-year duration of the simulation; and (ii) a five-year period followed by no further intervention, to simulate a single period of investment in population management. In each of these simulations, we modelled the interventions as (i) continuous (i.e. a constant rate of e.g. neutering) and (ii) annual (i.e. intervention applied once per year). Interventions were run at low, medium, and high coverages ( Table 2). As the processes are not equivalent, we apply different percentages for the intervention coverage (culling/neutering/sheltering) and the percent increase/decrease in parameter rates for the responsible ownership intervention. Intervention coverage refers to the proportion of dogs that are culled/neutered/sheltered per year (i.e. 20%, 40% and 70% annually) and, for responsible ownership interventions, the decrease in abandonment/roaming rate and increase in the adoption rate of shelter dogs (30%, 60% and 90% increase/decrease from baseline values). To model a low (20%), medium (40%) and high (70%) proportion of free-roaming dogs caught, but where half of the dogs were sheltered and half were neuteredand-returned, combined CNR and sheltering interventions were simulated at half-coverage (e.g. intervention rate of 0.7 was simulated by 0.35 neutered and 0.35 sheltered). For continuous interventions, sheltering (σ), culling (χ), and CNR (φ) were applied continuously during the length of the intervention. For annual interventions, σ, χ, and φ were applied to the ordinary differential equations using a forcing function applied at 12-month intervals. In simulations that included responsible ownership interventions, the decrease in owned dog abandonment/roaming (α) and the increase in shelter adoption (β) was assumed instantaneous and continuous (i.e. rates did not change throughout the intervention).
Model outputs. The primary outcome of interest was the impact of interventions on free-roaming dog population size. For interventions applied for the duration of the simulation, we calculated: (i) equilibrium population size for each population; (ii) percent decrease in free-roaming dog population; (iii) costs of intervention in terms of staff-time; and (iv) an overall welfare score. For interventions applied for a five-year period, we also calculated: (v) minimum free-roaming dog population size and percent reduction from initial population size; and (vi) the length of time between the end of the intervention and time-point at which the free-roaming dog population reached above 20,000 dogs (the assumed initial free-roaming dog population size of Lviv, based on our previous research 32 , see Supplementary Information for detail). The costs of population management interventions vary by country (e.g. staff salaries vary between countries) and by the method of application (e.g. method of culling, or resources provided in a shelter). To enable a comparison of the resources required for each intervention, the staff time (staff working-months) required to achieve the intervention coverage was calculated. While this does not incorporate the full costs of an intervention, as equipment (e.g. surgical equipment), advertising campaigns, travel costs for the animal care team, and facilities (e.g. clinic or shelter costs) are not included, it can be used as a proxy for intervention cost. Using data provided from VIER PFOTEN International, we estimated the average number of staff required to catch and neuter the free-roaming dog population and to house the shelter dog population in each intervention, using this data as a proxy for catching and sheltering/culling. The number of dogs that can be cared for per shelter staff varies by shelter. To account for this, we estimated two staff-to-dog ratios (low and high). Table 3 describes the staff requirements for the different interventions. www.nature.com/scientificreports/ Using the projected population sizes, the staff time required for each staff type (e.g. number of veterinarianmonths of work required) was calculated for each intervention. Relative salaries for the different staff types were estimated ( Table 3). The relative salaries were used to calculate the cost of the interventions by: [Staff time required × relative salary ] × €20,000. Where €20,000 was the estimated annual salary of a European veterinarian, allowing relative staff-time costs to be compared between the different interventions. Average annual costs were reported.
To provide overall welfare scores for each of the interventions, we apply the estimated welfare scores on a one to five scale, for each of the dog subpopulations, as determined by Hogasen et al. (2013) 22 . This scale is based on the Five Freedoms (freedom from hunger and thirst; freedom from discomfort; freedom from pain, injury, or disease; freedom to express normal behaviour; freedom from fear and distress 42,43 ) and was calculated using expert opinions from 60 veterinarians in Italy 22 . The scores were weighted by the participants' self-reported knowledge of different dog subpopulations, which resulted in the following scores: 2.8 for shelter dogs (W H ); 3.5 for owned dogs (W O ); 3.1 for neutered free-roaming dogs (W N ); and 2.3 for intact free-roaming dogs (W I ) 22 .
Using these estimated welfare scores, we calculated an average welfare score for the total dog population based on the model's projected population sizes for each subpopulation (Eq. 10). For interventions running for the duration of the simulation, the welfare score was calculated at the time point (t) when the population reached an equilibrium size. For interventions running for five years, the welfare score was calculated at the end of the five-year intervention. The percentage change in welfare scores from the baseline simulation were reported.

Model validation and sensitivity analysis.
A global sensitivity analysis was conducted on all parameters described in the baseline simulation and all interventions applied continuously, at high coverage, for the full duration of the simulation. A Latin square design algorithm was used in package "FME" 44 to sample the parameters within their range of values (Table 1). For the global sensitivity analysis on interventions, all parameter values were varied, apart from the parameters involved in the intervention (e.g. culling, neutering, abandonment/roaming rates). The effects of altering individual parameters (local sensitivity analysis) on the population equilibrium was also examined for the baseline simulation using the Latin square design algorithm to sample each parameter, individually, within their range of values. Sensitivity analyses were run for 100 simulations over 50 years solved with 0.01 step sizes.

Results
Model validation and sensitivity analysis in baseline simulation. The baseline simulation resulted in equilibrium sizes of 23,651 free-roaming dogs, 2086 shelter dogs, and 98,358 owned dogs. The population sizes settle at an equilibrium close, but not equal to, the carrying capacities for the free-roaming or owned populations due to the flows of abandonment/roaming (α), relinquishment (γ), free-roaming adoption (δ), and shelter adoption (β). In the global sensitivity analysis, for the baseline simulation, the simulated populations stabilised at an equilibrium between (A) 14,267 and 35,011 free-roaming dogs; (B) 596 and 4140 shelter dogs; (C) 84,737 and 109,236 owned dogs (Fig. S2). Results of the local sensitivity analysis are presented in Figs. S3 and S4. Table 4 outlines the impact of interventions applied for the full duration of the simulation. All interventions reduced the freeroaming dog population size (Fig. 2). Combined CNR and responsible ownership had the greatest reduction in free-roaming dog equilibrium population size at high (annual = 89%, continuous = 90% reduction) and medium (annual = 55%, continuous = 57% reduction) coverages. Culling most greatly reduced the population at low coverages (annual = 29%, continuous = 31%), although both sheltering and combined CNR and responsible ownership had similar effects. At all coverages, responsible ownership and CNR applied alone had only small effects on the free-roaming dog population, whereas culling and sheltering both had large effects. Combined CNR and sheltering was more effective than CNR, but less effective than sheltering applied alone.

Impacts of interventions applied for full duration of simulation (70 years).
The equilibrium population size fluctuated between annually applied interventions. Annual interventions that targeted removal from the population (i.e. culling and sheltering) fluctuated more than those that targeted recruitment (i.e. responsible ownership and CNR). When comparing the continuously applied interventions with the average of the minimum and maximum equilibrium population size for annual interventions, the resulting www.nature.com/scientificreports/ equilibrium population sizes were almost equal. Overall, the greatest proportion of neutered dogs was 0.77, achieved by annually applied high coverage combined CNR and responsible ownership (Fig. 3).
When applied continuously at a high coverage, responsible ownership was the cheapest intervention (average annual costs €120,690 for high, and €761,314 for low staff-to-dog ratio; Table 4), but only slightly cheaper than combined CNR and responsible ownership (€130,727 for high, and €762,187 for low staff-to-dog ratio). Combined CNR and responsible ownership was cheaper than CNR applied alone. This is due to the effectiveness of the combined method at reducing free-roaming dog population size, resulting in fewer dogs being neutered to maintain the neutering coverage. Sheltering applied continuously at high coverages was most expensive (€849,490 for high, and €7,753,805 for low staff-to-dog ratio). For interventions applied annually (responsible ownership was not included as continuous effects were assumed), culling was the cheapest intervention for all coverages Table 4. Impact of interventions on free-roaming (S), shelter (H) and owned (O) dog population equilibrium levels and percent change from baseline population equilibrium levels for interventions applied for the duration of the simulation at low, medium, and high coverages. The intervention with the greatest reduction in free-roaming dog population size for each periodicity-coverage combination is highlighted in italics. The greatest increase in welfare scores and the lowest average annual costs for each periodicity-coverage combination are highlighted in bold. *As annual interventions did not reach a single stable equilibrium point but varied between two points due to the annual interventions, the average of the minimum and maximum equilibrium population size was reported. www.nature.com/scientificreports/ for high staff-to-dog ratios, and combined CNR and responsible ownership was the cheapest intervention for all coverages for low staff-to-dog ratios. The greatest improvement in overall welfare score was achieved by combined CNR and responsible ownership for interventions applied continuously and annually at all coverages (Table 4), improving the welfare score by 6.58% and 6.71% respectively. The lowest improvement in welfare score was by responsible ownership at low coverages, which improved the welfare score by 0.57%.   Table 5 outlines the results for interventions applied for five years. For interventions applied annually and continuously, culling achieved the greatest reduction in free-roaming dog population size at all coverages (low: annual 34%, continuous 29%; medium: annual 57%, continuous 51%; high: annual 77%, continuous 71%), though sheltering was almost equally effective (Fig. 2).

Results for interventions applied for a single five-year period.
The free-roaming dog population returned most slowly to above 20,000 dogs in combined CNR and responsible ownership interventions for all coverages and periodicities (Table 5). For example, at high continuous coverages, the free-roaming dog population reached above 20,000 dogs after 5.3 years, compared to 1.1, 3.7, and 3.9 years for responsible ownership, sheltering and culling respectively. The free-roaming dog population Table 5. Impact of five-year intervention on minimum free-roaming dog population size and time taken between end of intervention and free-roaming dog population size reaching 20,000 dogs. The intervention with the overall longest time to reach baseline population size is highlighted in italics. The intervention with the greatest reduction in free-roaming dog population size, the longest time to reach baseline population size, the greatest increase in welfare score, and the lowest cost for each periodicity-coverage combination are highlighted in bold. www.nature.com/scientificreports/ returned more slowly to above 20,000 dogs for interventions with high coverages (i.e. due to higher magnitudes of effect). There was little difference between annual and continuous interventions. The intervention with the lowest total cost over the 5-year intervention was annually applied combined CNR and responsible ownership at a low coverage (€556,958 at high and €1,094,419 at low staff-to-dog ratios; Table 5). The most expensive intervention was sheltering applied annually at a high coverage (€6,326,746 at high and €57,130,659 at low staff-to-dog ratios).
The greatest improvement in overall welfare scores for both annual and continuous interventions was achieved by combined CNR and responsible ownership when applied at a high coverage (increase of 5.90% for both annual and continuous; Table 5). Responsible ownership applied continuously had the least improvement in welfare score when applied at a low coverage. This reduced the welfare score by 0.53%.

Discussion
This study provides important information on the predicted effectiveness of dog population management interventions whilst (i) explicitly incorporating the interactions between the different subpopulations, (ii) modelling a range of management methods at different periodicities and coverages, and (iii) also providing estimates of the welfare impacts and costs. Our results indicate that methods targeting multiple sources of population increase, such as combined CNR and responsible ownership interventions, may be most effective at reducing free-roaming dog population size when applied over longer periods of time, whilst also remaining cost-effective and maintaining highest levels of overall dog welfare. Sheltering and culling may be more effective over shorter periods of time but, compared to methods including neutering, the population may rapidly return to an equilibrium population size once management has ended. This model provides a template for assessing the impact of free-roaming dog population management in multiple contexts (e.g. rural, semi-rural, or urban populations) by applying subpopulation-specific parameters.
Reducing free-roaming dog population size can reduce the risks associated with free-roaming dogs (e.g. public health, livestock production, wildlife conservation). The results of this study find that high coverage and long-term combined CNR and responsible ownership may most greatly reduce these risks. Combined CNR and responsible ownership had a synergistic effect: CNR and responsible ownership applied in isolation were least effective at reducing free-roaming dog population size. These findings are similar to those reported in previous modelling studies that also suggest a synergistic effect when combining fertility control with reduced owned dog abandonment 15,23 . The synergistic effect can be explained by this method reducing multiple flows that can replenish the population size; without targeting multiple flows, the potential for the population to increase in size remains. Targeting multiple flows by combining responsible ownership interventions with other methods, such as culling or sheltering, may also have had strong reductive effects on the free-roaming dog population size. For example, combined culling with responsible ownership would increase the free-roaming dog death rate and reduce the flow of abandoned/roaming owned dogs to the free-roaming dog population. Indeed, other combinations of methods may also have been effective, such as combined CNR and culling. However, it is important to consider the logistical challenges in combining different methods. For example, combined CNR and culling may be challenging to apply in practice, as neutered dogs would need to be recognised and not culled to ensure efforts were not wasted. In this study, we only investigated methods that would feasibly be applied and have often been reported 27 .
When considering the individually applied interventions, culling and sheltering were most effective at reducing the free-roaming dog population size at all intervention coverages and periodicities (Fig. 2). Both culling and sheltering increase the rate of removal from the free-roaming dog population; culling by increasing the death rate and sheltering by increasing the sheltering rate. CNR and responsible ownership both reduce recruitment to the free-roaming dog population; CNR by reducing birth rates and responsible ownership by reducing owned dog abandonment and roaming rates. Our results suggest that, when considering interventions applied in isolation, increasing removal rates are more effective at reducing free-roaming dog population sizes than reducing recruitment rates. Our findings are similar to those reported by Amaku et al. 17 , who also reported culling to be more effective than CNR applied alone, but contradict those reported by Yoak et al. 25 and Hogasen et al. 22 who predicted CNR to be more effective than culling or sheltering at reducing free-roaming dog population size. This difference may be due to the differences in the subpopulations that were modelled, and the mathematical modelling methods and assumptions. For example, Yoak et al. 25 only modelled free-roaming dog population dynamics, with little abandonment of owned dogs, as appropriate to the dog population dynamics in Jaipur. If rates of abandonment were higher in Jaipur, the effects of CNR in Yoak et al. 25 's study may have been dampened. Yoak et al. 25 also implemented an agent-based modelling approach, compared to the differential equations applied by Amaku et al. and in our study, which may have resulted in the differences in management effectiveness.
Although our results suggest that CNR applied alone would be one of the least effective methods, there has been increasing interest over the last 20 years in the use of CNR to manage free-roaming dog populations 27 . Empirical studies have reported CNR to be effective in reducing population sizes between 12% over 1.5 years and 40% over 12 years at varying reported intervention coverages 18,45,46 . Our results suggest CNR may reduce the free-roaming dog population size up to 36% (high coverage) over 5 years, though it is hard to compare these results to empirical studies due to a lack of reporting of intervention coverage and management length (see 27 for review). It is also important to note that most empirical studies have been conducted in parts of India 18,45,46 , where rates of movement between the subpopulations (e.g. abandonment from owned to free-roaming dog populations) may vary compared to the parameters estimated in this study for Lviv. In areas where birth rates of free-roaming dogs are higher, CNR applied alone may be more effective at reducing the free-roaming dog population size.
There has been little focus on the impact and effectiveness of responsible ownership campaigns on freeroaming dog population size or human behaviour change 27  www.nature.com/scientificreports/ the abandonment and roaming of owned dogs could be as important as directing efforts to reduce the birth rate or increase the removal rate in areas where owned dog abandonment is a dominant source of free-roaming dogs. These findings are in line with those in previous modelling studies reporting synergistic effects of combined fertility control and restricted movement of owned dogs 26 , and emphasising the potential impact of owned dog abandonment dampening the effectiveness of CNR interventions 15,17,22 . Future dog population management efforts should aim to identify the sources of population increase (i.e. abandonment, roaming of owned dogs, or immigration) to understand the potential impact of CNR and other interventions in specific subpopulations. For example, combined CNR and responsible ownership would be more effective in areas where owned dog abandonment and roaming are higher, but would be less effective in areas with low levels of owned dog abandonment and roaming. Despite the promising potential effects of combined CNR and responsible ownership, the evidence is lacking on whether responsible ownership interventions can actually decrease abandonment and roaming rates and increase shelter adoption rates to the degree explored in this study. For example, we modelled a 90% reduction in abandonment/roaming rates, and a 90% increase in shelter adoption rates for the high intervention coverage. This would require a substantial increase in responsible ownership practices, and it is unclear whether this would be achievable in practice. Our lack of understanding of the effect of responsible ownership campaigns may be due to challenges in quantifying dog ownership practices. It is particularly challenging to accurately quantify the rate of owned dog abandonment, as this is likely to be under-reported 34 . Future studies may consider quantifying abandonment and roaming rates through questionnaire surveys, focus groups, or using local shelter relinquishment figures as an indicator of abandonment. It is important to measure the impact of responsible ownership campaigns, in combination with CNR efforts, to identify the most effective strategies in reducing abandonment rates and increasing responsible ownership practices.
In simulations where interventions were applied for five years, culling was most effective at reducing the freeroaming dog population size at all management coverages. In the present study, both CNR and combined CNR and responsible ownership were more effective at maintaining the free-roaming dog population at a lower size for longer than other interventions in the five-year intervention simulation (i.e. it took longer for the population to reach above 20,000 dogs after management ended). Individuals that have been removed through culling or sheltering may quickly be replaced through births or immigration, resulting in the population quickly returning to the carrying capacity. The longer lasting effects of neutering may be explained by the neutered individuals remaining in the population and, as they do not contribute to birth, dampening the potential for population growth.
Unsurprisingly, interventions including sheltering were most expensive overall (Tables 4 and 5). Sheltering increases the shelter dog population size and, if the rehoming rates stay the same, greater staff resources are required to care for the shelter dog population. Killing of shelter dogs would reduce the shelter dog population and thus reduce costs, but we did not model this as several countries have a no kill policy for shelter dogs 47 , only allowing euthanasia for behavioural or health problems. In these countries, sheltering as a method of population control may result in an increase in the shelter dog population and, without an improvement in rehoming rates, has the potential to lead to life-long stays in shelters and overcrowding. As this method is costly, it is potentially only a feasible option for higher income countries 10 .
The cheapest intervention overall in terms of staff costs was responsible ownership campaigns. Responsible ownership campaigns do not require dog catchers, veterinarians, or veterinary nurses. Although when applied alone this method was the least effective, the combination of CNR and responsible ownership was more effective at reducing the population size and was only marginally more costly than responsible ownership alone. As this method was effective in reducing population size, fewer dog catchers, veterinarians and nurses were required throughout the simulation, making it more cost-efficient in the long-term. Culling and CNR were more similar in costs, in terms of staff resources, than might be anticipated (Tables 4 and 5). Whilst CNR required higherpaid workers (veterinarians and veterinary nurses), as the simulation progressed, fewer intact dogs needed to be neutered over time to maintain the neutering coverage. This meant that fewer veterinarians were required over time, reducing the overall costs. Whereas with culling, similar numbers of dogs needed to be culled each year to reach the intervention coverage. This required a higher number of dog catchers to maintain the intervention coverage throughout the simulation, resulting in there being less of a difference between culling and CNR in terms of staff resources needed. It is worth noting that, in this study, we assumed that no veterinarians were required to cull dogs. This may not be reflective of reality, as veterinarians may be required to oversee culling in some countries. If we had included the costs a veterinarians for culling in the analyses, culling would have been less cost-effective. The costs of culling estimated in this study may therefore be lower than required in reality.
Our study has limitations in the cost analysis, as this did not encompass the full cost of the interventions, as it excluded staff training (e.g. in catching dogs, performing surgeries), facilities (e.g. clinics, kennels), and equipment costs. For example, creating and buying advertising space to promote a responsible ownership campaign may cost tens of thousands of euros (VIER PFOTEN International, personal communication), and this is not incorporated into the analysis. The cost of management is rarely provided by studies investigating the impact of dog population management 27 . Where costs do exist, they exist in disparate datasets: differing by country (i.e. country specific costs) and methods of application. For example, for CNR, costs may be lower for interventions that involve less staff-training and shorter pre-and post-surgical holding times, but may compromise dog health and welfare as a result 48 . Similarly, costs of running a shelter will vary depending on the shelter environment and management: shelters with larger dog-pen sizes, less enrichment, and lower staff-to-dog ratios may be less costly, but may also compromise health and welfare 49 . There are currently no data available to quantify the costs required for a responsible ownership campaign to reduce the abandonment and roaming of owned dogs by a quantified amount. To fully assess and compare the costs of different interventions requires integration of multiscale datasets across the dog population management field. Whilst this study does not encompass the full  11 . The greatest improvement in welfare score, compared to the baseline simulation, was achieved by combined CNR and responsible ownership. This is due to an overall reduction in free-roaming dog population size and an increase in the proportion of free-roaming dogs that were neutered, whose welfare is rated more highly than intact free-roaming dogs. However, as we used aggregated welfare scores, weighted by self-reported knowledge of dog subpopulations provided by veterinarians in Italy 22 , the scores may not be applicable to all countries. For example, the welfare of the free-roaming dog population in one country may be greater than another, due to country-specific risks to dog health and welfare. It is challenging to compare the welfare impact due to a lack of comparable welfare data within and between subpopulations 27 . Using this overall welfare score allows us to compare the potential welfare impact of different interventions on overall welfare. Dog population management methods may also have a short-term impact on dog health and welfare, depending on the method used. For example, there are important welfare risks related to CNR interventions, depending on dog handling procedures, standards of surgery, and post-operative care 48,50 . The impact of culling on dog welfare also depends on the method used. Historically, culling methods have been inhumane 51,52 . In some locations, mass killing is carried out by distributing poisons, most commonly strychnine 51 which is not a recommended method of killing from an animal welfare perspective 10 . In locations where poisoning is not used, the methods to restrain the dogs may also be of welfare concern. For example, methods of physical restraint include the holding of body parts by rope or metal tongs which can cause laceration or tissue damage 51 . Once restrained, methods to kill dogs have included electrocution, carbon monoxide poisoning and drowning 51 . These methods are not recommended by the World Organisation for Animal Health (OIE) but are still employed in some areas 10 . These additional welfare concerns should be considered when determining appropriate dog population management intervention.
A primary motive for dog population management in many countries globally is to reduce the risk of rabies transmission from dogs to humans 53 . Dogs are the maintenance host of the rabies virus and account for 99% of human-rabies transmissions 53 . Rabies causes an estimated 60,000 deaths annually 54 and is a disease that disproportionately affects lower-income countries 55 . Though population management for infectious disease control was not the focus of this study, our model could be extended to include infectious disease dynamics to explore the impact of vaccination strategies combined with dog population management interventions. The different subpopulations may have differences in accessibility for vaccination and the movement of dogs between the different subpopulations could have implications for disease control strategies. Incorporating dog subpopulation dynamics in infectious disease modelling could therefore be beneficial for identifying effective and efficient disease control strategies.
Several assumptions were made by this model. Firstly, the model did not differentiate between different age and sex categories. All individuals were assumed to contribute equally to the dynamics occurring between the different subpopulations. Currently the data is lacking at the level necessary to model the flows between subpopulations for different age and sex categories. Further study to determine these rates would be beneficial for informing future understanding of the dynamics between dog subpopulations. The model also did not include the effects of immigration and emigration explicitly. If interventions occur in neighbourhoods, instead of homogeneously throughout a city, there is the potential that this could increase movement of dogs from other parts of the city. We assumed that methods could be applied at consistent levels of coverage, regardless of the free-roaming dog population size. In reality, it may be more challenging to find dogs when the free-roaming dog population size is smaller.
The results of this study primarily apply to free-roaming dog populations in Europe. Populations in other geographic locations may differ in their movement within and between the subpopulations. Cultural factors within countries, such as religious beliefs, may affect public attitudes and acceptance of management methods. In areas with lower abandonment rates and/or opposition to neutering, combined CNR and responsible ownership interventions may be less effective. Despite this, this model highlights the importance of identifying sources of population increase and can be used by applying population specific parameters to identify effective and efficient management methods.

Conclusions
System dynamics modelling is a useful tool for investigating the potential impact of different dog population management methods. Overall, combined CNR and responsible ownership may be the most effective at reducing free-roaming dog population size, maintaining the population at a lower size, whilst also improving the overall welfare score of dogs, at low costs in terms of staff resources. Future dog population management would benefit from identifying and targeting all potential flows that may increase the free-roaming dog population size (such as abandonment, roaming owned dogs, and immigration). More studies are required to measure the impact of interventions, such as responsible ownership campaigns, on human behaviour change.